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We explore various properties of classical one-dimensional Wigner solids in the presence of disorder 
at T = in the context of a recently discovered Anderson transition of plasma modes in the random 
potential system. The extent to which the Wigner lattice is really a "crystal" rather than an 
amorphous solid is discussed for two types of disorder. We investigate the way in which zero point 
quantum fluctuations that would normally destroy the long range positional order are affected by 
the disorder. The probability density of particle spacings is examined analytically within a weak 
disorder approximation and compared to numerical calculations for two different realizations of 
disorder. We also discuss the multifractal character of the eigenmodes, the compressibility of the 
electrons, and the AC conductivity. 
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I. INTRODUCTION 

In 1958, E. Wigner predicted that an electron gas in 
a positive jellium background would, at low density, so- 
lidify into an ordered lattice^. Moreover, in classical sys- 
tems that have been studied in certain soft condensed 
matter models, the strong Coulomb repulsion is sufficient 
to cause a system of particles to crystallize^. 

In such a phase, the most natural physical quantities 
to consider are the elementary excitations. The charged 
collective modes or plasmons at T = are entities of this 
sort, the essential features of which are captured by a 
classical descriptioiii^. It is important to emphasize that 
the inclusion of strong disorder either chemically or, in 
the case of a low dimensional system, via a rough external 
substrate will result in the destruction of a well-defined 
base lattice. That is, the term "Wigner crystal" is no 
longer entirely descriptive in the presence of sufficiently 
strong disorder. Rather, that system is closer to a "glass" 
than a crystalline solid. 

Formally, the existence of Bragg peaks indicates that 
the system is a solid. However, others have discussed 
the possibility that the system retains quasi-long range 
order as a consequence of disorder—. Subtle implications 
are entailed by the term "glass" in describing a system in 
terms of its magnetic and topological properties. Here, 
we restrict our definition of "glass-like" to a system with 
a random positional configuration at zero magnetic field, 
even though a ID system cannot truly be topologically 
disordered as in higher dimensional systems^. 

In this article we focus on two aspects of one di- 
mensional disordered Wigner crystal(lDWC) systems — 
in particular the plasma oscillations and the statistical 
behavior of the equilibrium electron locations. We nu- 
merically study the plasma modes of the disordered sys- 
tem and calculate the compressibility and the AC con- 
ductivity. Regarding the structural properties, we con- 
sider whether the distribution of the spacings between ad- 
jacent electrons can be analytically understood in terms 
of the randomness in other parameters in the Hamilto- 
nian. A mathematical formalism is developed that pro- 



vides some insight into the connection between power 
spectrum of the randomness and the positional statis- 
tics. We also attempt to address the long-standing issue 
of whether a classical Wigner crystal in the presence of 
either random field or random charge disorder preserves 
some degree of crystalline integrity. 



II. THE STATISTICAL CONFIGURATION OF 
THE PARTICLES UNDER THE INFLUENCE OF 
WEAK DISORDER 

Let us now consider how disorder can affect the po- 
sitional ordering of the particles in the IDWC system. 
Evidently, the nature of disorder in such a system is not 
unique. It can exist in the spacings between the particle 
sites, the masses, the effective interaction between par- 
ticles, or one can have a set of charges pinned by a ran- 
dom potential. It is also the case that the introduction 
of randomness in one parameter can induce randomness 
in another. For example, randomness in particle-particle 
interactions or an external field may well result in the 
destruction of a periodic base lattice at mechanical equi- 
librium. 

The ordered system of size L is defined by the Hamil- 
tonian. 
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(2.1) 



where the quantities pi and m* refer to the momen- 
tum and effective mass of the i^^ particle, respectively. 
Our system lives in the low density regime in which the 
Coloumb interactions are much more influential than the 
kinetic energy, resulting in the crystalline ordering of the 
particles with the spatial coordinates {xi} separated with 
lattice constant a. Moreover, the system is stabilized by 
a positive Jellium neutralizing background. We wish to 
examine the vibrational modes of the system taking into 
account the full electrostatic force acting on each charge. 

One can add to this Hamiltonian a term representing 
the interaction of the charges with a random electrostatic 
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potential V{x), defined as, 

N 

V{x) — J a„ cos(27mx) + 6„ sin(27rrix) (2.2) 

n 

where the parameter J is a dimensionless coupling con- 
stant. The random variables a„ and 6„ have a mean equal 
to zero. Moreover, their distributions are controlled by 
independent Gaussian distributions the widths of which 
are independendent of n, yielding a "white noise" power 
spectrum. We will call this system Model A. 

Additionally, we can explicitly inject randomness into 
the single charge values Qi . We term this random charge 
system Model B, where the Hamiltonian is obtained from 
(|2.1|) by replacing the Cg with Qi x Qj. Some possible 
realizations of this model may be found in soft condensed 
matter systems such as in biopolymer arrays that have 
contain randomly assigned chemical constituents and un- 
screened Coulomb interactions^. 

In both Model A and Model B, the charges in our 
lattice will relax to mechanical equilibrium, at locations 
that differ from the evenly-spaced positions they occupy 
in the ordered case. Here, we characterize key properties 
of the equilibrium configuration, focusing on the statis- 
tical characteristics of the spacings between neighboring 
points. Finally, we discuss the extent to which the re- 
sulting lattice exhibits crystalline order. 

As a first step in answering these questions for Model 
B, we examine the product density functions for dif- 
ferent distributions of charges. This is shown in Ap- 
pendix [XI There is a direct relationship between this 
distribution and the distribution of forces in the case of 
strongly-screened interactions. The general distribution 
of the forces on a particular charge when interactions 
are unscreened requires a more involved analysis. Given 
this probability distribution, we turn to the distribution 
function of the spacings between the charges at equilib- 
rium. It turns out that an analytical result that is valid 
in the weak disorder regime — based on the assumption 
that the width of the single charge distributions is very 
small — agrees remarkably well with the results of numer- 
ical calculations. We will call this result the "viscosity 
approximation" in that it is based on the notion that 
the charges will deviate from the ordered lattice by dis- 
tances that are proportional to the total force acting on 
each charge cx Axi when they are equally spaced. 
That is, P{Fl"*') oc P{Axi), as if the system were placed 
in a highly viscous medium and initial displacements of 
the charges from their initial positions were predictive 
of their eventual positions. Given this assumption, one 
can work out the solution to Model B for all three of 
the distributions noted above starting with the gaussian 



system. 

The distribution of forces on the i^^ charge is given by 
the convolution: 



Pf{F,) = J PE{E,)Pq{qi)5{Eq - F,)dE,dq, 
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(2.3) 



We must first determine PsiEi), of which describes the 
likelihood that a given charge will have an electric field 
acting on it. The electric field acting on the i*"^ charge is 
Ptotai = E where ~ {xjx^y^ - Subsequently, if 

we examine the first charge in the array, the distribution 
of electric fields acting on it is a distribution governed by 
the convolution: 

OO OO OO 

PEiEi)= J dq2 J dqs--- J dqNPqiq2)PM---PqiqN) 

— OO — OO — OO 

X (^(926.1 + 936.1 H 1- qN^N.i - £^1) 

(2.4) 

This convolution can be simplified by making use of the 
Fourier representation of the Dirac delta function, 

^(926,1 + 936a H 1- <1N^N,1 - El) 

00 

^ ^ / exp[zfl(g26,i + 936,1 H ^ Qn^na - Ei)]dQ 

—00 

(2.5) 

More generally, for any given charge, the distribution of 
the electric fields PsiEi) is transformed into the following 
single integration. 
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where, we have integrated out the charges by performing 
an integral that defines g{k), 
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The resulting distribution of electric fields is in closed 
form, 
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The distribution of the force can now be determined by 
using ()2.3p . As a consequence of the periodic boundary 



conditions we can take, E^^i ^ks = Moreover, we can 
apply the exact sum E^^i l/'^^ = in the hmit of a 
large system size to the £,1 ^ summation, reducing ([2~ 
to a Gaussian distribution centered at zero: 



exp (-3 {E,f I (ttV 
(7r3a/3)'/' 



(2.9) 



The resulting force distribution is simply a convolution 
of two Gaussian distributions, 



exp (-3 {F.lq^'' I {-kg) - {q, - q^f /2a) 



(^V/3)^/'|ft| 

(2.10) 

We now turn to the other distribution functions. The 
system with uniform random charges, defined by the den- 
sity function (jA3p can also be evaluated. We can apply 
convolution of Eq. (12. 3p . Similar steps lead to the anal- 
ysis of Eq. ^^ri\ . 

f exp[if^gfcCfc,j] , 
^^-(""^ = J W,-W, 

Wi 

^ 2exp[z»a,,(W^2 + Wi)/2] X sin[n£k4W2 ~ Wi)/2] 

{W2-Wl)nCk,^ 

After applying eauation (|2.4p we have: 

00 

PEiE.,) = l- f e-^"^'x 

ZTT 



(2.11) 



e»»a-(^^+^^)/^sin[»a,,(M^2 - l^^i)/2] | 



n£k4W2-Wi)/2 

Similarily, for binary disorder we have 

00 



(2.12) 



(2.13) 



This leads to the final force distributions for the uni- 
form and the binary systems respectively as: 



priF,) 



Ei i-iriFjWi) - Ei {~-inFjW2) 
27r(VF2 - Wi) 



^ e»"^''-W+^^)/^sin[»a,,(Vt/2 - T4^i)/2] ^ 



n£kAW2-Wi)/2 



(2.14) 



and 



dq^ pt,in(^^) ^ 



/ exp(-^^]j^,/Ql) expj^inFjQ^] 
V 47rQi 



TV 

X I []e^^^«''-(Qi+Q^V2cos[l]a..:(Q2-gi)/2] | dl^ 

yk^i 

(2.15) 

where Ei{x) is the exponential integral defined as: 



Ei{x) 
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(2.16) 



The density of particle spacings is simply proportional to 
these force distributions. 

This derivation in terms of probability convolutions 
can similarly be extended to the random potential sys- 
tem of Model A. We previously defined the random po- 
tential in Eq. (|2.2p . To simplify the analysis we will 
only consider the cosine contributions and will perform 
the analysis for a Gaussian distribution of the random 
amplitudes. The resulting Fourier decomposition yields 
the expression, 

1 ^ 

V{x) = —= ^ a„ cos(n7ra;/L) (2.17) 
The force on the i"^ charge is then given by, 

L 

F[xi) = -cqV' {x)\,j.^ = —q^ a„(»7r) sin(n7ra:,/L) 



n=l 



(2.18) 

It is evident that the source of the randomness is not 
a coupling constant between particles as was the case in 
Model B. Rather, we relate probability density of random 
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Fourier amplitudes, a„ ("white noise"), to the particle ity convolution for the force is given as: 
probability density of successive spacings. The probabil- 
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P{Fi) = eo j dai j da2 ■ ■ ■ j daLPa{ai)Pa{a2) 
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(2.19) 



■Pa{aL)5 ( jjj:^ ( ai sin(-^) + 2a2 sin(— ^) . . . Lol sin(7r.Ti) - Fi 



r 



As a consequence of Eq. (|2.19p having the same math- 
ematical structure as (|2.4p . one can carry out the same 
steps that lead to Eq. I|2.8p . yielding the expression: 



eoexp (-Ffl (2a ^ti (^sin(l^; 



2(tX]^=i (x37^sin(n7rxi/i) 



(2.20) 

Our final result for the distribution of nearest neigh- 
bor spacings of Model A is simply a Gaussian distribution 
with a width that depends on the relative spatial loca- 
tion, Xi of any given particle. This corresponds to the 
same Gaussian distribution type as the random ampli- 
tudes an used to construct the random potential. How- 
ever, the width differs, as the spacings distribution con- 
tains a width that is a function of the precise character- 
istics of the random potential. 

In Model B, one can naively compare all of the derived 
probability densities of nearest neighbor spacings to the 
probability densities of force couplings given earlier by 
Eqs. (HHni), (Emi), (EUSl). The binary system contains 
three characteristic peaks in both quantities, showing 
similar behavior. Furthermore, in the Gaussian case both 
the force couplings and the nearest neighbor spacings fol- 
low a convolution of two Gaussian distributions. On the 
other hand, for the uniform distribution the probability 
density of nearest neighbor spacings does not follow the 
logarithmic behavior of the force couplings. Ultimately, 
this indicates the importance of the distribution type in 
detailing the equilibrium structural properties of IDWC 
systems. 

As remarked earlier, the viscosity approximation as- 
sumes that the particles will deviate from the ordered 
lattice by a displacement that is proportional to the to- 
tal force on each particle when the particles are equally 
spaced. This approximation seems reasonable when the 
disorder is relatively weak compared to the interaction 
strength of the Coloumb force. However, it is necessary 
test this hypothesis by performing a numerical relaxation 
of the particles for both Models A and B in finite sized 
systems. We detail the procedure and compare the nu- 
merical results to our previously derived expressions in 



the next section. 



III. THE NUMERICAL RELAXATION OF THE 
PARTICLES TO EQUILIBRIUM 

The numerical relaxation of a disordered IDWC to its 
minimum energy configuration requires the development 
of two major components. First one must address the 
complications that arise from performing computations 
on finite size systems with long ranged interactions. Ev- 
idently, the surface effects, as a result of long range in- 
teractions in finite sized systems are not negligible. One 
possible route to dealing with this complication is to en- 
force periodic boundary conditions. The standard ap- 
proach involves the Ewald summation technique^. We 
have derived our version of the Ewald potential displayed 
as Eq. (|B2p in Appendix [B] that represents the inter- 
action energy between charges on a periodic image, of 
which is infinitely repeated to assert the periodic bound- 
ary conditions. However, the resulting rapidly converg- 
ing sums still demand computationally intensive efforts. 
We have discovered that the summation can be repro- 
duced in closed form in terms of a simple function that 
reproduces the correct short distance behavior and has 
the correct periodicity. Consider the lattice summation 
of the \/R Couloumb potential representing the infinite 
sum of n images of length d. The interaction potential 
between two particles separated by a distance x within a 
image is given by the expression 



E 



(-1)" 



(n + x/d) 



TT |csc(7rx/d)| (3.1) 



The result is a simple trigonometric function, of which is 
easy to implement numerically. FiglT] shows a graphical 
comparison of the expressions (|B2[) and (|3.ip . 

Next we consider the other required numerical tool. 
An energy minimization algorithm is neccessary to find 
the spatial coordinates at which the particles arrange 
themselves given a particular realization of disorder in 
the Hamiltonian. These coordinates satisfy a system of 
linear equations, namely the simultaneous zeros of the 
total force on each charge. A very powerful technique is 
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FIG. 1: A comparison of the exact and Ewald representations 
of the Coulomb interaction energy 



FIG. 3: A comparison of the analytical and numerical distri- 
bution of nearest beighbor spacings -(continued) 
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FIG. 2: A comparison of the analytical and numerical distri- 
bution of nearest neighbor spacings 



the Newton-Raphson(NR) method^. In the NR method 
one recursively computes the Hessian or second deriva- 
tive matrix and multiplies its inverse by the initial po- 
sitions. The resulting product is used to increment the 
coordinates towards the equilibrium configuration. One 
interesting property of the Hessian matrix is that it is 
precisely the Dynamical Matrix one uses to compute the 
phonon spectrum of a vibrating lattice. In our IDWC 
system, we can examine the plasma oscillations by simply 
diagonalizing the Hessian of a relaxed array of particles. 
Using the interaction potential of Eq. (|3.ip , the Hessian 
of a finite size lattice has the forn*^ 



D(R R') = Sn,R' 



R" 



£c=R-R" 



dx^ 



a;=R-R' 

(3.2) 

This yields an N x N matrix that can one use to carry 
out the standard phonon analysis. We have implemented 
numerical relaxations on chains of length N = 1024, for 
each of the random charge and random field distributions 
discussed in the previous section. Numerical uncertainty 
exists in the amount of residual total force acting on a 
particular charge. Therefore, the charges were relaxed 
so that the residual forces were only about « 10~^ in 



relative magnitude (scaled to unity), at which point our 
system is quite close to equilibrium. We have compared 
the theoretical predictions of the viscosity approximation 
given by equations (PH)) . (^13)) . ^JU^i, and (P?^ to the 
numerical data in Figs. [2l[3l We find excellent agreement 
between the viscosity approximation and numerical cal- 
culations in terms in the general behavior of the resulting 
distributions. Furthermore, it is clear that the statistical 
arrangement of the disorder has an explicit dependence 
on the distribution type and width. The correlation is 
quite striking in the case of the binary random charge sys- 
tem, where the three characteristic peaks are mirrored in 
both the numerical and analytical plots. This suggests 
that an experimental probe of the structural arrange- 
ment of disordered systems with long-ranged interactions 
would reveal a pronounced difference between the chemi- 
cally disordered binary system and the random potential 
system, of which represents a form of substrate disorder. 
It is also worth noting that the effect could also exist in 
several soft condensed matter models in which Wigner 
crystal ordering has been discussed in the context of bi- 
ological systems^. 



IV. CRYSTALLINE ORDER? 

It is well known that there is a strict absence of long 
range crystalline order for systems with short-ranged 
interactions that exist in dimensions d < 2, resulting 
from the prevailing influence of thermal and quantum 
fluctuations. This notion has been explored by vari- 
ous authors^. For Wigner crystals in particular, the un- 
screened Coulomb interactions factor into these melting 
criteria through the precise wavevector dependence of the 
longitudinal and transverse vibrational modes that lead 
to diverging values of the Lindemann parameter {u^} in 
a low enough dimensionality. 

As we have demonstrated in the preceding section, 
the inclusion of quenched disorder induces randomness 
in particle spacings. Furthermore, another main conse- 
quence of the disorder is that eigenmodes are not plane 
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waves. However, as we will show in our investigation 
of the Bragg peaks, the Fourier space is not necessarily 
destroyed. 

In our IDWC system, let us isolate the role of quenched 
disorder in making the particle array more amorphous 
while neglecting the known effects of stochastic zero point 
disturbances studied by other authors^^. Again, this is 
valid in systems in which Wigner-like ordering is present 
but those entropic contributions are suppressed. A more 
useful route to the exploration of the possibility of long 
range order involves the structure factor 5^, which can 
be defined in terms of the equilibrium positions of the ith 
particle, a;^"^": 



1,2,3- •• (4.1) 



The crystalline character of the system will manifest itself 
through the existence of at least one delta function, or 
Bragg peak, in the spectrum of Sk within some particular 
scale in Fourier space. 

From the exact dispersion relatior>ii*i^ of the clean 
ID Wigner crystal, uj{k) ^ |fc| log^^^(l//c), it has been 
shown that the correlation functions decay much slower 

than any power law, (\u{x) — u(0)]^^ oc ■y/ln(a;). The 

resulting prediction of strong Break peaks in a scattering 
experiment follows from the scattering intensity tending 

like, I{x) oc Exp[-A([u{x) ~ u{Q)f'^], of which is rel- 
atively slow compared to the typical, strong exponen- 
tial decay observed in systems with short-ranged inter- 
actions. According to the literature, systems possessing 
this sort of decay are usually labeled as having "quasi- 
long range order"—. 

The quantity Sk was calculated for numerically relaxed 
ensembles in both models A and B at different relative 
disorder strengths. In the case of Model A, the disorder 
strength strength is the dimensionless combination / J 
for values of J defined in Eq. (|2.2p ). while for model B, 
it is simply the single charge distribution width W = 
1 + Wo{W2 — Wi), where Wq = is the ordered case. 

The Sk data is shown in Figs. [l]and[5l The plots sug- 
gest a difference between both models in terms of their 
structural properties. It appears that both models re- 
tain their crystalline composition for a regime of disor- 
der strengths. Moreover, in the high disorder regime, the 
Bragg peaks are lost. The details of the transition will 
be examined more carefully in a future work. A notable 
difference between the systems is that in model A the 
Bragg peaks strongly resemble delta functions while in 
model B, the scattering intensity contains an exponen- 
tially broadened shoulder tending on both sides of the 
peak at higher disorder strengths, although the presence 
of peak shoulders in Model B is not evident at low dis- 
order. Hence for weak disorder both models are crystals 
based on the strictest definition requiring the existence 
of just one Bragg peak. 

In order to appropriately make sense of Figs. 2] and 
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FIG. 4: Sk for different relative interaction strengths, model 
A, L—256. At high k, or low disorder strength, in the delocal- 
ization critical regime, the system maintains the composition 
of the first Bragg peak, indicating the preservation of some 
sort of crystalline aspect to the system. 
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FIG. 5: Sk for different disorder strengths, model B, L=256. 



[S] one should focus on the first Bragg peak and how 
it behaves as the disorder strength is increased. Ev- 
idently, the peak height decreases for both models as 
function of disorder strength. It is also worth noting 
that the peaks decrease across Fourier space accord- 
ing to a Debye- Waller damping envelope, tending like, 
Idw ^ exp (— ((g7i(0))^)). The precise reduction of the 
first Bragg peak height g was studied for both models 
and we determined a power law behavior of the peak 
height in terms of the disorder strength. We performed a 
linear fit shown in Fig. [5] on a logarithmic scale yield- 



for values of 



ing qaIk) ^ K"" and gsiWo) 
a « -1.7272 ± 0.0696 and /3 « -0.9217 ± 0.0324. Ap- 
parently, this power law behavior breaks down at higher 
disorder strengths, as the system becomes more amor- 
phous. 
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FIG. 6: The behavior of the height of the first Bragg peak g 
for different disorder strengths. 



FIG. 7: Averaged participation moments showing multifractal 
scahng, Model A 



V. PLASMA OSCILLATIONS OF THE 
DISORDERED SYSTEM 

The localization properties of the plasma eigenstates 
of models A and B have been explored in detail by the 
authors of this paper—, where they reported a vibra- 
tional delocalization transition at lower eigenfrequencies 
in Model A, while no transition observed in Model B. 
Additionally, they performed a finite size scaling analysis 
that yielded a scaling behavior of the relative interaction 
strength / J and an approximate value of the corre- 
lation length exponent. In this section we will focus on 
other aspects of the plasma eigenstates. In particular, 
we will explore multifractal properties of the participa- 
tion moments, the electronic compressibility of the pure 
and disordered systems, and the frequency dependent be- 
havior of the AC conductivity. 

An important physical quantity associated with spa- 
tially localized eigenfunctions is the participation ratio 
P^— of a given moment 5, defined as: 



q 



1 



k=l 



2q 



(5.1) 



where w^. is the amplitude of the jth plasma eigenmode 
at site k. 

For this particular quantity one must account for the 
numerical and statistical uncertainties on finite length 
chains that scale as 1/VL. One can smooth out the 
fiuctuations of P-j by summing over all the eigenstates 
and ensemble-averaging over several realizations of dis- 
order (represented by the over line bar). This yields the 
following length dependent quantity for finite-sized sys- 
tems: 



PqiL) = - 



1 



p3 



(5.2) 



For the 3D A nderson Model the critical scaling prop- 
erties of Pq{L) were derived by Wegner— using renor- 
malization group arguments in a 2 -I- e expansion of the 



2 particle spectral function of a system of interacting Q 
matrices. The common view is that all models exhibit- 
ing an Anderson transition will possess critical eigenfunc- 
tions that can be used to construct the quantity Pq{L), 
having a length dependence that scales like: 



PqiL) cx L^'f^^ 



(5.3) 



or in other words, Pq{L) exhibits multifractal scaling, in 
that one obtains a different exponent Dq for each mo- 
ment q. Other authors such as Lima et. al. have re- 
lated this scaling behavior to the wave-packet dynam- 
ics of the power-law bond disordered one-dimensional 
Anderson model with hopping amplitudes decreasing as 
Hnm (X |n — to| . They were able to extract a time- 
dependent scaling analysis of the participation moments 
by performing a finite size scaling analysis of the elec- 
tronic return probabilitj^. Furthermore, they report an 
asymptotic value of Doo = 0.5. We have performed a 
similar calculation shown in Fig. [T) Indeed, the critical 
eigenfunctions of Model A exhibit multifractal behavior. 
However the asymptotic value D^o w 0.2, differs from the 
value 0.5 obtained in the analogous electronic tightbind- 
ing system. 

It is important to note that the behavior of the plasma 
oscillations have important consequences for physical 
quantities such as the electron compressibility k. We re- 
mark that our definition of the k only includes the elec- 
tron gas itself and does not consider the compensating 
background. The main reason for using this definition 
is that the background contribution does not enter into 
the experimentally determined values of n when using 
the double well field penetration technique, of which has 
been widely applied to measure the compressibility of in- 
teracting electron systems in semiconductor heterostruc- 
ture systemsi^. The thermodynamic compressibility can 
be related to the volume fluctuations of a system around 
its average value, given by the relation 



K (X 



[Avy 



{V} 



(5.4) 
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One can associate the normal modes u{x) with these vol- 
ume fluctuations in a finite system by the relation, 



1 / A 1 ( du 




(5.5) 



where oj^ is the nth eigenfrequency and a;„ labels the nth 
oscillator. 

For the ordered — or clean — Wigner crystal, the precise 
behavior of the compressibility can be extracted by ex- 
ploiting the spatial periodicity of the longitudinal plasma 
eigenfrequencies ujl &t long wavelengthaii. 



hm „ 



(5.6) 



Apparently, for the classical ordered system one can sub- 
stitute the long wavelength form of the dispersion rela- 
tions uj{k) ~ \k\ log^^^(l/A;), into (|5.6|) to yield a vanish- 
ing compressibility, k = O^i. 

Similarly, the quantum calculations having the disper- 
sion relationi^ 



UJRPA{k) 



2eo 



(5.7) 



also yield k — 0, where f/ is the Fermi velocity and d is 
a finite system width. 

We also note that our definition of the compressibility 
does not always yield a vanishing result in the thermo- 
dynamic limit. For a harmonic oscillator system with 
random couplings and nearest neighbor interactions we 
observe a finite compressibility. This appears to be con- 
sistent with a dispersion Lo{k) that is linear in fc, in the 
ordered system, of which, k = via Eq. (|5.6p . This leads 
us to believe that the compressibility of the charges is de- 
pendent on the range of the interactions, not the nature 
of the disorder. 

We have also applied Eq. (|5.5p to Models A and B 
at various system sizes and have observed a vanishing 
compressibility k = 0, in the limit of a large system, 
for both cases. This is a significant result, given that 
the properties of disordered elastic systems has recently 
been subject of intense researchi^. We are confident that 
the methods used in this article can easily be extended 
to two dimensional systems, where a major unresolved a 
problem is the compressibility behavior observed in 2D SI 
MOSFETS. Various research teams have experimentally 
observed a vanishing compressibility on the insulating 
side of the well known 2D Metal-Insulator transition's^. 

There has been much speculation as to whether the 
insulating phase is an incompressible zero field Wigner 
crystal. Certain authors have suggested the use of the 
gaussian variational replica method(GSM) as a viable 
tool for understanding disordered Wigner crystals. How- 
ever, it is known that these elastic models assume a con- 
tinuum limit that washes out the full microscopic struc- 
tural arrangement of the electrons; for example one could 
not discuss the distribution of nearest neighbor spacings 
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FIG. 8: Re[(T((x;)] vs. uj for a fixed relative disorder strength 
in both models k — 0.1, Wo — 0.5 at different system sizes. 



in a random charge system like Model B using the GSM, 
given that the spacings between nearest neighbor parti- 
cles is quantity observed only by a discrete description as 
we have considered in this paper. Moreover, even if one 
wishes to apply a continuum description, one needs to 
determine the elastic constants to effectively utilize the 
GSM method for computing the relevant physical quanti- 
ties, thus rendering previous replica based studies incon- 
clusive. On the other hand, we have been able to calcu- 
late the compressibility by circumventing these field the- 
oretical limitations on elastic systems(constrained by the 
Larkin scale) through the use of simple classical phonon 
methods. We hope to discuss the 2D calculations in a 
future article, given the fact that we have also developed 
some precise mathematical tools for handling the 2D long 
range Coulomb interactions, of which are analogous to 
the exact summation in Eq. (j3.ip . 

Another important quantity that derives from the 
plasma oscillations is the AC conductivity. The plasmon 
propagator G{uj) is itself the direct response to an ex- 
ternal electric field interacting with the electrons. G{uj) 
can be determined from the Hessian or dynamical ma- 
trix, defined by Eq. I\3.'2\i , as constructed in the following 
resolvent: 



R R' 



ujH-D{R,R') 



(5.8) 



The AC conductivity is then given by the quantity^: 



a{Lj) = iAijjG{Lo) 



(5.9) 



The effects of disorder on the zero temperature AC 
transport properties was discussed by Giamarchi and 
Schultz^^ for ID electron systems and Chitra et al. for 2D 
Wigner Crystals^"^. With reference to this earlier work 
we have calculated the AC conductivity via Eq. (|5.9p . 
The finite size effects are shown in Fig|8] for both Mod- 
els A and B. The key quantities of interest are the peak 
height P and its associated frequency ojp. For the pure 
crystalline system corresponds to the Drude driving 
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frequency. The disorder acts as a broadening mechanism, 
by which Up is shifted to a new frequency known as the 
pinning frequency and also the peak height P is shifted 
as well. We have qualitatively studied how ujp scales with 
the relative interaction strength of the Coloumb potential 
over the random potential previously defined as k in the 
critical Model A. With reference to our earlier work, the 
observed interaction scaling in Model A can be utilized to 
relate the critical eigenfrequency at which delocalization 
takes place ujc to the pinning frequency Up. We have 
plotted the behavior of Re[cr(a;)] at various interaction 
strengths shown in Fig. [TOl There is a noticeable increase 
in peak height with increasing interaction strength and a 
subtle decrease in LUp as displayed in the window of Fig. 
1101 Consequently, this implies that for critical frequen- 
cies uJc closer to the upper band edge there is a smaller 
pinning frequency. This can be further interpreted as 
connecting stronger plasma eigenmode localization with 
higher pinning frequencies. We remark that more care- 
ful scaling studies, both analytical and numerical must 
performed for a deeper understanding of the connection 
between critical depinning transitions and vibrational de- 
localization. 

Another important property of the AC conductivity 
is the high frequency behavior of Re[o'(a;)]. We exam- 
ined the limit lu >> ujp and observed an lu dependence 
going like Re[(T(cj)] ^ l/oj^ for both Models A and B. 
This is shown in Fig. [51 Certain authors have deter- 
mined a universal lo^ dependence at low frequencies for 
ID disordered Wigner Crystals while others have estab- 
lished a for high frequencies from renormalization 
group studies of interacting ID electron systems, where 
K is an effective scaling coupling constant for repulsive 
interactions^. These earlier studies only consider short 
ranged interactions while our system takes into account 
the true long-ranged nature of the Coulomb interactions. 
The latter form would apply to our classical systems for 
the value K 1, yielding a dependence. For the 
purposes of comparison with similar ID disordered elec- 
tron models, our dependence on uj differs, for our results 
are precisely the same as the 2D system studied by Gia- 
marchi et al. using the GSM technique, where they also 
observe Re[fT(w)] ~ behavior for high frequencies. 



VI. CONCLUSION AND EXPERIMENTAL 
PROSPECTS 

In this paper we have considered a wide range of phe- 
nomena in ID disordered Wigner solids that follow from 
the behavior of the plasma oscillations. We have de- 
veloped a new mathematical formalism for studying the 
spatial arrangement of the particles at equilibrium as 
a consequence of a certain type of randomness in the 
Hamiltonian. We have applied probability convolutions 
to quantify the distribution of nearest neighbor particle 
spacings by assuming the particle coordinates will devi- 
ate from crystalline configuration by distances that are 
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FIG. 9; Lo^Re[a{uj)] vs. lu for a fixed relative disorder strength 
in both models k — 0.1, Wo = 0.5 at L=512. The zero slope 
regime at higher frequencies clearly suggests a l/u)^ depen- 
dence. 
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FIG. 10: Re[(T(tj)] vs. u at different relative disorder strength 
values, k(L=512). The inset shows the shift of the pinning 
frequency ujp as the disorder is decreased. 



proportional to the total force on each particle when the 
system is ordered. If one adiabatically "turns on" the 
randomness, one can imagine the charges propagating in 
a viscous medium, and if the random contribution to the 
forces on the charges is small a linear relaxation of the 
particles might be expected to hold as a good approxima- 
tion. We have tested this viscosity approximation against 
numerically relaxed ensembles, and observe what appears 
to be remarkably good agreement. A stronger statistical 
analysis will be performed in future papers that will fur- 
ther develop this probability convolution formalism on 
more rigorous grounds. 

We have also explored the basic question of whether 
a IDWC in the presence of quenched disorder retains 
its crystalline composition. We have argued that both 
a random potential system and a random charge system 
generate an internal stiffness that resists the system from 
completely melting in the thermodynamic limit, at least 
when the effect of randomness is not too great. 

We briefly remark that a melting transition has been 
observed at a critical value of disorder strength. However, 
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this will be formally reported in a future article. 

In the process of studying Models A and B numeri- 
cally, we have revised existing computational paradigms 
for incorporating long ranged Coulomb forces in systems 
with periodic boundary conditions and disorder. A nu- 
merical relaxation scheme has been outlined along with 
a procedure for studying the collective modes of disor- 
dered Wigner crystals or other "glasslike" structures in 
general. 

We conclude by emphasizing the necessity for experi- 
mental probes of the density response functions in order 
that one may investigate the possibility of a vibrational 
delocalization transition and our predictions of a zero 
compressibility phase. Glasson et. al. have already ob- 
served Wigner Crystal ordering in a quasi-lD system.^^. 
We consider the possibility of well-controlled manipu- 
lation of the substrate randomness as a viable avenue 
by which the scaling behavior of models A and B can 
be investigated experimentally. Furthermore, based on 
the observed connection between the delocalization tran- 
sition and the pinning frequency ujp, AC conductivity 
measurements are also quite essential for elucidating our 
understanding of the basic physics involved. 



After determining the range of values for t, we can 
solve for the joint distribution for a product of charges: 



(A4) 



As a more germane alternative, we also consider a binary 
distribution for Model B, defined as simply a sum of two 
delta functions with equal statistical weights. This par- 
ticular type of disorder can be interpreted as a form of 
chemical disorder or a type of binary alloy. We consider 
the binding of chemical constituents to the Wigner crys- 
tal such that there are effectively two values of charges 
random distributed throughout the system: 



1 



1 



P.iQ) ^ T^SiQ - Qi) + -S{Q - Q2) 



(A5) 



The product distribution simply consists of three peaks 
at the values Ql, Qi x Q2, and Q^- 

P{t) = -^5{Ql -t) + ^S{QiQ2 ^t) + ^S{Ql - t) (A6) 
Finally, we consider the gaussian product distribution: 
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APPENDIX A: THE DISTRIBUTION OF 
RANDOM CHARGE PRODUCTS 

In general, the charge product density function is given 
as a convolution of two different densities^^: 



Pit) = J j Pq4Q^)Pcl,,{QJ)5{Q,QJ-t)dQdQJ (Al) 

If the two charges share a common density function, then: 

t .dQ 



Pit) = / P,{Q)P,i 



Q'\Q\ 



(A2) 



The range of values for t must be carefully determined 
for various regimes, depending on the type and width of 
distribution Pq{Qi). We start by examining the case of a 
uniform distribution of random single charges in Model 
B, in which 



pm 







Wi<Q<W2 
otherwise 



(A3) 



Evidently, nuclcons and electrons can only have integer 
values of the charge e, so this particular distribution is 
unphysical and unrealistic. We consider this distribution 
only as a version of a randomly-distributed coupling con- 
stant and for the purposes of determining which features 
of the system are robust with respect to distribution type. 



P{t) 



iQiQj~t)dQidQj 
(A7) 



For fi = this integral can be evaluated exactly to yield 
P{t) = (^)i^o f 4) (A8) 

where Kq is a modified Bessel function of the second 
kind.— For a non-zero mean the product distribution can 
be determined numerically. 



APPENDIX B: THE EWALD SUMMATION OF A 
QUASI-ID PERIODIC SYSTEM 

We carry out the Ewald summation technique for a 
quasi- ID system of charges in a uniform neutralizing 
background. The electrostatic potential of interest has 
the following form: 



1 



dx/ L 



\r-nL\ J ^(r - xf + b'^ 
The first term can be expanded out as an integral: 



(Bl) 



\r — nL\ 







This integral can be broken up into two parts by intro- 
ducing a suitable cutoff ^: 
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The 2nd integral can be evaluated through a variable 
substitution and can be left as is: 



TTX 



oo roo 



E 
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Next, we apply the Poisson sum formula to the 1st inte- 
gral. 



Lastly we subtract the divergent contribution at m = 
arising from the neutralizing background, 



^(r-a;)2 + fe2 J J L 



If we include the other m = term from the previous 
summand the contributions to the sum become, 



E 



m— — oo 



g2™rmg-t(r-nL)^-l/2^^^^ 



The integral on n can carried out by a standard comple- 
tion of the squares trick yielding. 



E 4 



t 



It is important that we properly subtract the divergent 
contribution of the m = term. The resulting summand 
contains a lower incomplete gamma function: 



,dt 

T J ' t 



where C is some constant. Our final expression for the 
long-ranged electrostatic potential becomes, 



oo 

oo „ 

V{r)= 



n— — oo 



(B2) 



+ Ee'™^^r(0,^)-ln(;3) 



It is important to note that the expression converges to 
a value that is independent of the cutoff (3. 
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